Corona_Analysis

The 2019-2020 Coronavirus Pandemic Analysis

Contact: Smith Research

BACKGROUND & APPROACH

I wanted to track and trend the coronavirus outbreak on my own curiosity. There are some interesting questions that may fall out of this, as it is a very historic moment, including scientifically and analytically (we have a large amount of data being shared across the globe, analyzed in real-time). The world has come to a halt because of it.
This analysis attempts to answer the following questions (more to come):

1. What does the trend of the pandemic look like to date?
   
2. What are future case predictions based on historical model?
3. What interesting quirks or patterns emerge?

ASSUMPTIONS & LIMITATIONS: * This data is limited by the source. I realized early on that depending on source there were conflicting # of cases. Originally I was using JHU data… but this was always ‘ahead’ of the Our World In Data. I noticed that JHU’s website was buggy- you clicked on the U.S. stats but it didn’t reflect the U.S.. So I changed data sources to be more consistent with what is presented in the media (and Our World In Data has more extensive plots I can compare my own to). An interesting aside might be why the discrepancy? Was I missing something?
* Defintiions are important as is the idea that multiple varibales accumulate in things like total cases (more testing for example).

SOURCE RAW DATA: * https://ourworldindata.org/coronavirus
* https://github.com/CSSEGISandData/COVID-19/
*

INPUT DATA LOCATION: github (https://github.com/sbs87/coronavirus/tree/master/data)

OUTPUT DATA LOCATIOn: github (https://github.com/sbs87/coronavirus/tree/master/results)

TIMESTAMP

Start: ##—— Sun May 17 19:30:32 2020 ——##

PRE-ANALYSIS

The following sections are outside the scope of the ‘analysis’ but are still needed to prepare everything

UPSTREAM PROCESSING/ANALYSIS

1. Google Mobility Scraping, script available at get_google_mobility.py

# Mobility data has to be extracted from Google PDF reports using a web scraping script (python , written by Peter Simone, https://github.com/petersim1/MIT_COVID19)

# See get_google_mobility.py for local script 

python3 get_google_mobility.py
# writes csv file of mobility data as "mobility.csv"

SET UP ENVIORNMENT

Load libraries and set global variables

# timestamp start
timestamp()
## ##------ Sun May 17 19:30:32 2020 ------##

# clear previous enviornment
rm(list = ls())

##------------------------------------------
## LIBRARIES
##------------------------------------------
library(plyr)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  3.0.0     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(ggplot2)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(plot.utils)
library(utils)
library(knitr)

##------------------------------------------

##------------------------------------------
# GLOBAL VARIABLES
##------------------------------------------
user_name <- Sys.info()["user"]
working_dir <- paste0("/Users/", user_name, "/Projects/coronavirus/")  # don't forget trailing /
results_dir <- paste0(working_dir, "results/")  # assumes diretory exists
results_dir_custom <- paste0(results_dir, "custom/")  # assumes diretory exists


Corona_Cases.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
Corona_Cases.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv"
Corona_Deaths.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv"
Corona_Deaths.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"

Corona_Cases.fn <- paste0(working_dir, "data/", basename(Corona_Cases.source_url))
Corona_Cases.US.fn <- paste0(working_dir, "data/", basename(Corona_Cases.US.source_url))
Corona_Deaths.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.source_url))
Corona_Deaths.US.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.US.source_url))
default_theme <- theme_bw() + theme(text = element_text(size = 14))  # fix this
##------------------------------------------

FUNCTIONS

List of functions

function_name	description
prediction_model	outputs case estumate for given log-linear moder parameters slope and intercept
make_long	converts input data to long format (specialized cases)
name_overlaps	outputs the column names intersection and set diffs of two data frame
find_linear_index	finds the first date at which linearaity occurs

##------------------------------------------
## FUNCTION: prediction_model
##------------------------------------------
## --- //// ----
# Takes days vs log10 (case) linear model parameters and a set of days since 100 cases and outputs a dataframe with total number of predicted cases for those days
## --- //// ----
prediction_model<-function(m=1,b=0,days=1){
  total_cases<-m*days+b
  total_cases.log<-log(total_cases,10)
  prediction<-data.frame(days=days,Total_confirmed_cases_perstate=total_cases)
  return(prediction)
}
##------------------------------------------

##------------------------------------------
## FUNCTION: make_long
##------------------------------------------
## --- //// ----
# Takes wide-format case data and converts into long format, using date and total cases as variable/values. Also enforces standardization/assumes data struture naming by using fixed variable name, value name, id.vars, 
## --- //// ----
make_long<-function(data_in,variable.name = "Date",
                   value.name = "Total_confirmed_cases",
                   id.vars=c("case_type","Province.State","Country.Region","Lat","Long","City","Population")){

long_data<-melt(data_in,
                id.vars = id.vars,
                variable.name=variable.name,
                value.name=value.name)
return(long_data)

}
##------------------------------------------

## THIS WILL BE IN UTILS AT SOME POINT
name_overlaps<-function(df1,df2){
i<-intersect(names(df1),
names(df2))
sd1<-setdiff(names(df1),
names(df2))
sd2<-setdiff(names(df2),names(df1))
cat("intersection:\n",paste(i,"\n"))
cat("in df1 but not df2:\n",paste(sd1,"\n"))
cat("in df2 but not df1:\n",paste(sd2,"\n"))
return(list("int"=i,"sd_1_2"=sd1,"sd_2_1"=sd2))
}

##------------------------------------------

##------------------------------------------
## FUNCTION: find_linear_index
##------------------------------------------
## --- //// ----
# Find date at which total case data is linear (for a given data frame) 
## --- //// ----

find_linear_index<-function(tmp,running_avg=5){
  tmp$Total_confirmed_cases_perstate.log<-log(tmp$Total_confirmed_cases_perstate,2)
  derivative<-data.frame(matrix(nrow = nrow(tmp),ncol = 4))
  names(derivative)<-c("m.time","mm.time","cumsum","date")
  
  # First derivative
  for(t in 2:nrow(tmp)){
    slope.t<- tmp[t,"Total_confirmed_cases_perstate.log"]- tmp[t-1,"Total_confirmed_cases_perstate.log"]
    derivative[t,"m.time"]<-slope.t
    derivative[t,"date"]<-as.Date(tmp[t,"Date"])
  }
  
  # Second derivative
  for(t in 2:nrow(derivative)){
    slope.t<- derivative[t,"m.time"]- derivative[t-1,"m.time"]
    derivative[t,"mm.time"]<-slope.t
  }
  
  #Compute running sum of second derivative (window = 5). Choose point at which within 0.2
  for(t in running_avg:nrow(derivative)){
    slope.t<- sum(abs(derivative[t:(t-4),"mm.time"])<0.2,na.rm = T)
    derivative[t,"cumsum"]<-slope.t
  }
  
  #Find date -5 from the stablility point
  linear_begin<-min(derivative[!is.na(derivative$cumsum) & derivative$cumsum==running_avg,"date"])-running_avg
  
  return(linear_begin)
}

READ IN DATA

• total number of cases. current source: https://github.com/CSSEGISandData (precvious source https://ourworldindata.org/coronavirus)

# Q: do we want to archive previous versions? Maybe an auto git mv?

##------------------------------------------
## Download and read in latest data from github
##------------------------------------------
download.file(Corona_Cases.source_url, destfile = Corona_Cases.fn)
Corona_Totals.raw <- read.csv(Corona_Cases.fn, header = T, stringsAsFactors = F)

download.file(Corona_Cases.US.source_url, destfile = Corona_Cases.US.fn)
Corona_Totals.US.raw <- read.csv(Corona_Cases.US.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.source_url, destfile = Corona_Deaths.fn)
Corona_Deaths.raw <- read.csv(Corona_Deaths.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.US.source_url, destfile = Corona_Deaths.US.fn)
Corona_Deaths.US.raw <- read.csv(Corona_Deaths.US.fn, header = T, stringsAsFactors = F)

# latest date on all data:
paste("US deaths:", names(Corona_Deaths.US.raw)[ncol(Corona_Deaths.US.raw)])

## [1] "US deaths: X5.16.20"

paste("US total:", names(Corona_Totals.US.raw)[ncol(Corona_Totals.US.raw)])

## [1] "US total: X5.16.20"

paste("World deaths:", names(Corona_Deaths.raw)[ncol(Corona_Deaths.raw)])

## [1] "World deaths: X5.16.20"

paste("World total:", names(Corona_Totals.raw)[ncol(Corona_Totals.raw)])

## [1] "World total: X5.16.20"

PROCESS DATA

• Convert to long format
  
• Fix date formatting/convert to numeric date
  
• Log10 transform total # cases

##------------------------------------------
## Combine death and total data frames
##------------------------------------------
Corona_Totals.raw$case_type<-"total"
Corona_Totals.US.raw$case_type<-"total"
Corona_Deaths.raw$case_type<-"death"
Corona_Deaths.US.raw$case_type<-"death"

# for some reason, Population listed in US death file but not for other data... Weird. When combining, all datasets will have this column, but US deaths is the only useful one.  
Corona_Totals.US.raw$Population<-"NA" 
Corona_Totals.raw$Population<-"NA"
Corona_Deaths.raw$Population<-"NA"

Corona_Cases.raw<-rbind(Corona_Totals.raw,Corona_Deaths.raw)
Corona_Cases.US.raw<-rbind(Corona_Totals.US.raw,Corona_Deaths.US.raw)
#TODO: custom utils- setdiff, intersect names... option to output in merging too
##------------------------------------------
# prepare raw datasets for eventual combining
##------------------------------------------
Corona_Cases.raw$City<-"NA" # US-level data has Cities
Corona_Cases.US.raw$Country_Region<-"US_state" # To differentiate from World-level stats

Corona_Cases.US.raw<-plyr::rename(Corona_Cases.US.raw,c("Province_State"="Province.State",
                                                  "Country_Region"="Country.Region",
                                                  "Long_"="Long",
                                                  "Admin2"="City"))


##------------------------------------------
## Convert to long format
##------------------------------------------
#JHU has a gross file format. It's in wide format with each column is the date in MM/DD/YY. So read this in as raw data but trasnform it to be better suited for analysis
# Furthermore, the World and US level data is formatted differently, containing different columns, etc. Recitfy this and combine the world-level stats with U.S. level stats.

Corona_Cases.long<-rbind(make_long(select(Corona_Cases.US.raw,-c(UID,iso2,iso3,code3,FIPS,Combined_Key))),
make_long(Corona_Cases.raw))


##------------------------------------------
## Fix date formatting, convert to numeric date
##------------------------------------------
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "^X",replacement = "0") # leading 0 read in as X
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "20$",replacement = "2020") # ends in .20 and not 2020
Corona_Cases.long$Date<-as.Date(Corona_Cases.long$Date,format = "%m.%d.%y")
Corona_Cases.long$Date.numeric<-as.numeric(Corona_Cases.long$Date)

kable(table(select(Corona_Cases.long,c("Country.Region","case_type"))),caption = "Number of death and total case longitudinal datapoints per geographical region")

Number of death and total case longitudinal datapoints per geographical region

	death	total
Afghanistan	116	116
Albania	116	116
Algeria	116	116
Andorra	116	116
Angola	116	116
Antigua and Barbuda	116	116
Argentina	116	116
Armenia	116	116
Australia	928	928
Austria	116	116
Azerbaijan	116	116
Bahamas	116	116
Bahrain	116	116
Bangladesh	116	116
Barbados	116	116
Belarus	116	116
Belgium	116	116
Belize	116	116
Benin	116	116
Bhutan	116	116
Bolivia	116	116
Bosnia and Herzegovina	116	116
Botswana	116	116
Brazil	116	116
Brunei	116	116
Bulgaria	116	116
Burkina Faso	116	116
Burma	116	116
Burundi	116	116
Cabo Verde	116	116
Cambodia	116	116
Cameroon	116	116
Canada	1624	1624
Central African Republic	116	116
Chad	116	116
Chile	116	116
China	3828	3828
Colombia	116	116
Comoros	116	116
Congo (Brazzaville)	116	116
Congo (Kinshasa)	116	116
Costa Rica	116	116
Cote d’Ivoire	116	116
Croatia	116	116
Cuba	116	116
Cyprus	116	116
Czechia	116	116
Denmark	348	348
Diamond Princess	116	116
Djibouti	116	116
Dominica	116	116
Dominican Republic	116	116
Ecuador	116	116
Egypt	116	116
El Salvador	116	116
Equatorial Guinea	116	116
Eritrea	116	116
Estonia	116	116
Eswatini	116	116
Ethiopia	116	116
Fiji	116	116
Finland	116	116
France	1276	1276
Gabon	116	116
Gambia	116	116
Georgia	116	116
Germany	116	116
Ghana	116	116
Greece	116	116
Grenada	116	116
Guatemala	116	116
Guinea	116	116
Guinea-Bissau	116	116
Guyana	116	116
Haiti	116	116
Holy See	116	116
Honduras	116	116
Hungary	116	116
Iceland	116	116
India	116	116
Indonesia	116	116
Iran	116	116
Iraq	116	116
Ireland	116	116
Israel	116	116
Italy	116	116
Jamaica	116	116
Japan	116	116
Jordan	116	116
Kazakhstan	116	116
Kenya	116	116
Korea, South	116	116
Kosovo	116	116
Kuwait	116	116
Kyrgyzstan	116	116
Laos	116	116
Latvia	116	116
Lebanon	116	116
Lesotho	116	116
Liberia	116	116
Libya	116	116
Liechtenstein	116	116
Lithuania	116	116
Luxembourg	116	116
Madagascar	116	116
Malawi	116	116
Malaysia	116	116
Maldives	116	116
Mali	116	116
Malta	116	116
Mauritania	116	116
Mauritius	116	116
Mexico	116	116
Moldova	116	116
Monaco	116	116
Mongolia	116	116
Montenegro	116	116
Morocco	116	116
Mozambique	116	116
MS Zaandam	116	116
Namibia	116	116
Nepal	116	116
Netherlands	580	580
New Zealand	116	116
Nicaragua	116	116
Niger	116	116
Nigeria	116	116
North Macedonia	116	116
Norway	116	116
Oman	116	116
Pakistan	116	116
Panama	116	116
Papua New Guinea	116	116
Paraguay	116	116
Peru	116	116
Philippines	116	116
Poland	116	116
Portugal	116	116
Qatar	116	116
Romania	116	116
Russia	116	116
Rwanda	116	116
Saint Kitts and Nevis	116	116
Saint Lucia	116	116
Saint Vincent and the Grenadines	116	116
San Marino	116	116
Sao Tome and Principe	116	116
Saudi Arabia	116	116
Senegal	116	116
Serbia	116	116
Seychelles	116	116
Sierra Leone	116	116
Singapore	116	116
Slovakia	116	116
Slovenia	116	116
Somalia	116	116
South Africa	116	116
South Sudan	116	116
Spain	116	116
Sri Lanka	116	116
Sudan	116	116
Suriname	116	116
Sweden	116	116
Switzerland	116	116
Syria	116	116
Taiwan*	116	116
Tajikistan	116	116
Tanzania	116	116
Thailand	116	116
Timor-Leste	116	116
Togo	116	116
Trinidad and Tobago	116	116
Tunisia	116	116
Turkey	116	116
Uganda	116	116
Ukraine	116	116
United Arab Emirates	116	116
United Kingdom	1276	1276
Uruguay	116	116
US	116	116
US_state	378276	378276
Uzbekistan	116	116
Venezuela	116	116
Vietnam	116	116
West Bank and Gaza	116	116
Western Sahara	116	116
Yemen	116	116
Zambia	116	116
Zimbabwe	116	116

# Decouple population and lat/long data, refactor to make it more tidy
metadata_columns<-c("Lat","Long","Population")
metadata<-unique(select(filter(Corona_Cases.long,case_type=="death"),c("Country.Region","Province.State","City",all_of(metadata_columns))))
Corona_Cases.long<-select(Corona_Cases.long,-all_of(metadata_columns))

# Some counties are not summarized on the country level. collapse all but US
Corona_Cases.long<-rbind.fill(ddply(filter(Corona_Cases.long,!Country.Region=="US_state"),c("case_type","Country.Region","Date","Date.numeric"),summarise,Total_confirmed_cases=sum(Total_confirmed_cases)),filter(Corona_Cases.long,Country.Region=="US_state"))

# Put total case and deaths side-by-side (wide)
Corona_Cases<-spread(Corona_Cases.long,key = case_type,value = Total_confirmed_cases)

#Compute moratlity rate
Corona_Cases$mortality_rate<-Corona_Cases$death/Corona_Cases$total

#TMP
Corona_Cases<-plyr::rename(Corona_Cases,c("total"="Total_confirmed_cases","death"="Total_confirmed_deaths"))

##------------------------------------------
## log10 transform total # cases
##------------------------------------------
Corona_Cases$Total_confirmed_cases.log<-log(Corona_Cases$Total_confirmed_cases,10)
Corona_Cases$Total_confirmed_deaths.log<-log(Corona_Cases$Total_confirmed_deaths,10)
##------------------------------------------
       
##------------------------------------------
## Compute # of days since 100th for US data
##------------------------------------------

# Find day that 100th case was found for Country/Province. NOTE: Non US countries may have weird provinces. For example, Frane is summairzed at the country level but also had 3 providences. I've only ensured the U.S. case100 works... so the case100_date for U.S. is summarized both for the entire country (regardless of state) and on a per-state level. 
# TODO: consider city-level summary as well. This data may be sparse

Corona_Cases<-merge(Corona_Cases,ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,case100_date=min(Date.numeric)))
Corona_Cases$Days_since_100<-Corona_Cases$Date.numeric-Corona_Cases$case100_date

##------------------------------------------
## Add population and lat/long data (CURRENTLY US ONLY)
##------------------------------------------

kable(filter(metadata,(is.na(Country.Region) | is.na(Population) )) %>% select(c("Country.Region","Province.State","City")) %>% unique(),caption = "Regions for which either population or Country is NA")

Regions for which either population or Country is NA

Country.Region	Province.State	City

# Drop missing data 
metadata<-filter(metadata,!(is.na(Country.Region) | is.na(Population) ))
# Convert remaining pop to numeric
metadata$Population<-as.numeric(metadata$Population)

## Warning: NAs introduced by coercion

# Add metadata to cases
Corona_Cases<-merge(Corona_Cases,metadata,all.x = T)

##------------------------------------------
## Compute total and death cases relative to population 
##------------------------------------------

Corona_Cases$Total_confirmed_cases.per100<-100*Corona_Cases$Total_confirmed_cases/Corona_Cases$Population
Corona_Cases$Total_confirmed_deaths.per100<-100*Corona_Cases$Total_confirmed_deaths/Corona_Cases$Population


##------------------------------------------
## Filter df for US state-wide stats
##------------------------------------------

Corona_Cases.US_state<-filter(Corona_Cases,Country.Region=="US_state" & Total_confirmed_cases>0 ) 
kable(table(select(Corona_Cases.US_state,c("Province.State"))),caption = "Number of longitudinal datapoints (total/death) per state")

Number of longitudinal datapoints (total/death) per state

Var1	Freq
Alabama	3596
Alaska	642
Arizona	926
Arkansas	3705
California	3545
Colorado	3249
Connecticut	547
Delaware	223
Diamond Princess	61
District of Columbia	62
Florida	3886
Georgia	8561
Grand Princess	62
Guam	62
Hawaii	322
Idaho	1694
Illinois	4788
Indiana	4997
Iowa	4393
Kansas	3573
Kentucky	5299
Louisiana	3630
Maine	917
Maryland	1433
Massachusetts	945
Michigan	4334
Minnesota	3999
Mississippi	4495
Missouri	4831
Montana	1520
Nebraska	2628
Nevada	667
New Hampshire	642
New Jersey	1406
New Mexico	1461
New York	3511
North Carolina	5238
North Dakota	1679
Northern Mariana Islands	47
Ohio	4685
Oklahoma	3468
Oregon	1811
Pennsylvania	3730
Puerto Rico	62
Rhode Island	364
South Carolina	2639
South Dakota	2166
Tennessee	5043
Texas	10292
Utah	863
Vermont	855
Virgin Islands	62
Virginia	6552
Washington	2410
West Virginia	2321
Wisconsin	3490
Wyoming	1081

Corona_Cases.US_state<-merge(Corona_Cases.US_state,ddply(filter(Corona_Cases.US_state,Total_confirmed_cases>100),c("Province.State"),summarise,case100_date_state=min(Date.numeric)))
Corona_Cases.US_state$Days_since_100_state<-Corona_Cases.US_state$Date.numeric-Corona_Cases.US_state$case100_date_state

ANALYSIS

Q1: What is the trend in cases, mortality across geopgraphical regions?

Plot # of cases vs time
* For each geographical set:
* comparative longitudinal case trend (absolute & log scale)
* comparative longitudinal mortality trend
* death vs total correlation

question	dataset	x	y	color	facet	pch	dimentions
comparative_longitudinal_case_trend	long	time	log_cases	geography	none (case type?)	case_type	[15, 50, 4] geography x (2 scale?) case type
comparative longitudinal case trend	long	time	cases	geography	case_type	?	[15, 50, 4] geography x (2+ scale) case type
comparative longitudinal mortality trend	wide	time	mortality rate	geography	none	none	[15, 50, 4] geography
death vs total correlation	wide	cases	deaths	geography	none	none	[15, 50, 4] geography

# total cases vs time
# death cases vs time
# mortality rate vs time
# death vs mortality


  # death vs mortality
  # total & death case vs time (same plot)

#<question> <x> <y> <colored> <facet> <dataset>
## trend in case/deaths over time, comapred across regions <time> <log cases> <geography*> <none> <.wide>
## trend in case/deaths over time, comapred across regions <time> <cases> <geography*> <case_type> <.long>
## trend in mortality rate over time, comapred across regions <time> <mortality rate> <geography*> <none>
## how are death/mortality related/correlated? <time> <log cases> <geography*> <none>
## how are death and case load correlated? <cases> <deaths>

# lm for each?? - > apply lm from each region starting from 100th case. m, b associated with each.
    # input: geographical regsion, logcase vs day (100th case)
    # output: m, b for each geographical region ID



#total/death on same plot-  diffeer by 2 logs, so when plotting log, use pch. when plotting absolute, need to use free scales
#when plotting death and case on same, melt. 

#CoronaCases - > filter sets (3)
  #world - choose countries with sufficent data

N<-ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,n=length(Country.Region))
ggplot(filter(N,n<100),aes(x=n))+
  geom_histogram()+
  default_theme+
  ggtitle("Distribution of number of days with at least 100 confirmed cases for each region")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

kable(arrange(N,-n),caption="Sorted number of days with at least 100 confirmed cases")

Sorted number of days with at least 100 confirmed cases

Country.Region	n
US_state	28900
China	116
Diamond Princess	97
Korea, South	87
Japan	86
Italy	84
Iran	81
Singapore	78
France	77
Germany	77
Spain	76
US	75
Switzerland	73
United Kingdom	73
Belgium	72
Netherlands	72
Norway	72
Sweden	72
Austria	70
Malaysia	69
Australia	68
Bahrain	68
Denmark	68
Canada	67
Qatar	67
Iceland	66
Brazil	65
Czechia	65
Finland	65
Greece	65
Iraq	65
Israel	65
Portugal	65
Slovenia	65
Egypt	64
Estonia	64
India	64
Ireland	64
Kuwait	64
Philippines	64
Poland	64
Romania	64
Saudi Arabia	64
Indonesia	63
Lebanon	63
San Marino	63
Thailand	63
Chile	62
Pakistan	62
Luxembourg	61
Peru	61
Russia	61
Ecuador	60
Mexico	60
Slovakia	60
South Africa	60
United Arab Emirates	60
Armenia	59
Colombia	59
Croatia	59
Panama	59
Serbia	59
Taiwan*	59
Turkey	59
Argentina	58
Bulgaria	58
Latvia	58
Uruguay	58
Algeria	57
Costa Rica	57
Dominican Republic	57
Hungary	57
Andorra	56
Bosnia and Herzegovina	56
Jordan	56
Lithuania	56
Morocco	56
New Zealand	56
North Macedonia	56
Vietnam	56
Albania	55
Cyprus	55
Malta	55
Moldova	55
Brunei	54
Burkina Faso	54
Sri Lanka	54
Tunisia	54
Ukraine	53
Azerbaijan	52
Ghana	52
Kazakhstan	52
Oman	52
Senegal	52
Venezuela	52
Afghanistan	51
Cote d’Ivoire	51
Cuba	50
Mauritius	50
Uzbekistan	50
Cambodia	49
Cameroon	49
Honduras	49
Nigeria	49
West Bank and Gaza	49
Belarus	48
Georgia	48
Bolivia	47
Kosovo	47
Kyrgyzstan	47
Montenegro	47
Congo (Kinshasa)	46
Kenya	45
Niger	44
Guinea	43
Rwanda	43
Trinidad and Tobago	43
Paraguay	42
Bangladesh	41
Djibouti	39
El Salvador	38
Guatemala	37
Madagascar	36
Mali	35
Congo (Brazzaville)	32
Jamaica	32
Gabon	30
Somalia	30
Tanzania	30
Ethiopia	29
Burma	28
Sudan	27
Liberia	26
Maldives	24
Equatorial Guinea	23
Cabo Verde	21
Sierra Leone	19
Guinea-Bissau	18
Togo	18
Zambia	17
Eswatini	16
Chad	15
Tajikistan	14
Haiti	12
Sao Tome and Principe	12
Benin	10
Nepal	10
Uganda	10
Central African Republic	9
South Sudan	9
Guyana	7
Mozambique	6
Yemen	2
Mongolia	1

# Pick top 15 countries with data
max_colors<-12
# find way to fix this- China has diff provences. Plot doesnt look right...
sufficient_data<-arrange(filter(N,!Country.Region %in% c("US_state", "Diamond Princess")),-n)[1:max_colors,]
kable(sufficient_data,caption = paste0("Top ",max_colors," countries with sufficient data"))

Top 12 countries with sufficient data

Country.Region	n
China	116
Korea, South	87
Japan	86
Italy	84
Iran	81
Singapore	78
France	77
Germany	77
Spain	76
US	75
Switzerland	73
United Kingdom	73

Corona_Cases.world<-filter(Corona_Cases,Country.Region %in% c(sufficient_data$Country.Region))


  #us 
  #    - by state
Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
# summarize 
#!City %in% c("Unassigned") 
  #    - specific cities
#mortality_rate!=Inf & mortality_rate<=1
Corona_Cases.UScity<-filter(Corona_Cases,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") & City %in% c("Bucks","Baltimore City", "New York","Burlington"))

measure_vars_long<-c("Total_confirmed_cases.log","Total_confirmed_cases","Total_confirmed_deaths","Total_confirmed_deaths.log")
melt_arg_list<-list(variable.name = "case_type",value.name = "cases",measure.vars = c("Total_confirmed_cases","Total_confirmed_deaths"))
melt_arg_list$data=NULL


melt_arg_list$data=select(Corona_Cases.world,-ends_with(match = "log"))
Corona_Cases.world.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.UScity,-ends_with(match = "log"))
Corona_Cases.UScity.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.US_state,-ends_with(match = "log"))
Corona_Cases.US_state.long<-do.call(melt,melt_arg_list)

Corona_Cases.world.long$cases.log<-log(Corona_Cases.world.long$cases,10)
Corona_Cases.US_state.long$cases.log<-log(Corona_Cases.US_state.long$cases,10)
Corona_Cases.UScity.long$cases.log<-log(Corona_Cases.UScity.long$cases,10)


# what is the current death and total case load for US? For world? For states?
#-absolute
#-log

# what is mortality rate (US, world)
#-absolute

#how is death and case correlated? (US, world)
#-absolute

#Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
#Corona_Cases.US.case100<-filter(Corona_Cases.US, Days_since_100>=0)
# linear model parameters
#(model_fit<-lm(formula = Total_confirmed_cases.log~Days_since_100,data= Corona_Cases.US.case100 ))

#(slope<-model_fit$coefficients[2])
#(intercept<-model_fit$coefficients[1])

# Correlation coefficient
#cor(x = Corona_Cases.US.case100$Days_since_100,y = Corona_Cases.US.case100$Total_confirmed_cases.log)

##------------------------------------------
## Plot World Data
##------------------------------------------
# Timestamp for world
timestamp_plot.world<-paste("Most recent date for which data available:",max(Corona_Cases.world$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")


# Base template for plots
baseplot.world<-ggplot(data=NULL,aes(x=Days_since_100,col=Country.Region))+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))


##/////////////////////////
### Plot Longitudinal cases

(Corona_Cases.world.long.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world)
    )

(Corona_Cases.world.loglong.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases.log))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases.log))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world))

##/////////////////////////
### Plot Longitudinal mortality rate

(Corona_Cases.world.mortality.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world,aes(y=mortality_rate))+
    geom_line(data=Corona_Cases.world,aes(y=mortality_rate))+
    ylim(c(0,0.3))+
    ggtitle(timestamp_plot.world))

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.world.casecor.plot<-ggplot(Corona_Cases.world,aes(x=Total_confirmed_cases,y=Total_confirmed_deaths,col=Country.Region))+
  geom_point()+
  geom_line()+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
    ggtitle(timestamp_plot.world))

### Write polots

write_plot(Corona_Cases.world.long.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.long.plot.png"

write_plot(Corona_Cases.world.loglong.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.loglong.plot.png"

write_plot(Corona_Cases.world.mortality.plot,wd = results_dir)

## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.mortality.plot.png"

write_plot(Corona_Cases.world.casecor.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.casecor.plot.png"

##------------------------------------------
## Plot US State Data
##-----------------------------------------

baseplot.US<-ggplot(data=NULL,aes(x=Days_since_100_state,col=case_type))+
  default_theme+
  facet_wrap(~Province.State)+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))

Corona_Cases.US_state.long.plot<-baseplot.US+geom_point(data=Corona_Cases.US_state.long,aes(y=cases.log))
##------------------------------------------
## Plot US City Data
##-----------------------------------------

Corona_Cases.US.plotdata<-filter(Corona_Cases.US_state,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") &
                                   City %in% c("Bucks","Baltimore City", "New York","Burlington") &
                                   Total_confirmed_cases>0) 
timestamp_plot<-paste("Most recent date for which data available:",max(Corona_Cases.US.plotdata$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")

city_colors<-c("Bucks"='#beaed4',"Baltimore City"='#386cb0', "New York"='#7fc97f',"Burlington"='#fdc086')

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.city.loglong.plot<-ggplot(melt(Corona_Cases.US.plotdata,measure.vars = c("Total_confirmed_cases.log","Total_confirmed_deaths.log"),variable.name = "case_type",value.name = "cases"),aes(x=Date,y=cases,col=City,pch=case_type))+
  geom_point(size=4)+
    geom_line()+
  default_theme+
  #facet_wrap(~case_type)+
    ggtitle(paste("Log10 total and death cases over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
    scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.long.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State,scales = "free_y")+
  ggtitle(paste("MD, PA, NJ total cases over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))
+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.mortality.plot<-ggplot(Corona_Cases.US.plotdata,aes(x=Date,y=mortality_rate,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Mortality rate (deaths/total) over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.casecor.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(y=Total_confirmed_deaths,x=Total_confirmed_cases,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Correlation of death vs total cases,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
  scale_color_manual(values = city_colors))

(Corona_Cases.city.long.normalized.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases.per100,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State)+
  ggtitle(paste("MD, PA, NJ total cases over time per 100 people,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)  +
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

write_plot(Corona_Cases.city.long.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.plot.png"

write_plot(Corona_Cases.city.loglong.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.loglong.plot.png"

write_plot(Corona_Cases.city.mortality.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.mortality.plot.png"

write_plot(Corona_Cases.city.casecor.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.casecor.plot.png"

write_plot(Corona_Cases.city.long.normalized.plot,wd = results_dir_custom)

## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.normalized.plot.png"

Q1b what is the model

Fit the cases to a linear model 1. Find time at which the case vs date becomes linear in each plot
2. Fit linear model for each city

# What is the predict # of cases for the next few days?
# How is the model performing historically?

Corona_Cases.US_state.summary<-ddply(Corona_Cases.US_state,
                                     c("Province.State","Date"),
                                     summarise,
                                     Total_confirmed_cases_perstate=sum(Total_confirmed_cases)) %>% 
    filter(Total_confirmed_cases_perstate>100)

# Compute the states with the most cases (for coloring and for linear model)
top_states_totals<-head(ddply(Corona_Cases.US_state.summary,c("Province.State"),summarise, Total_confirmed_cases_perstate.max=max(Total_confirmed_cases_perstate)) %>% arrange(-Total_confirmed_cases_perstate.max),n=max_colors)

kable(top_states_totals,caption = "Top 12 States, total count ")
top_states<-top_states_totals$Province.State

# Manually fix states so that Maryland is switched out for New York
top_states_modified<-c(top_states[top_states !="New York"],"Maryland")

# Plot with all states:
(Corona_Cases.US_state.summary.plot<-ggplot(Corona_Cases.US_state.summary,aes(x=Date,y=Total_confirmed_cases_perstate))+
  geom_point()+
  geom_point(data=filter(Corona_Cases.US_state.summary,Province.State %in% top_states),aes(col=Province.State))+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  default_theme+
  theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
  ggtitle("Total confirmed cases per state, top 12 colored")+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Fit linear model to time vs total cases
##-----------------------------------------

# First, find the date at which each state's cases vs time becomes lienar (2nd derivative is about 0)
li<-ddply(Corona_Cases.US_state.summary,c("Province.State"),find_linear_index)

# Compute linear model for each state starting at the point at which data becomes linear
for(i in 1:nrow(li)){
  Province.State.i<-li[i,"Province.State"]
  date.i<-li[i,"V1"]
  data.i<-filter(Corona_Cases.US_state.summary,Province.State==Province.State.i & as.numeric(Date) >= date.i)
  model_results<-lm(data.i,formula = Total_confirmed_cases_perstate~Date)
  slope<-model_results$coefficients[2]
  intercept<-model_results$coefficients[1]
  li[li$Province.State==Province.State.i,"m"]<-slope
  li[li$Province.State==Province.State.i,"b"]<-intercept
  }

# Compute top state case load with fitted model

(Corona_Cases.US_state.lm.plot<-ggplot(filter(Corona_Cases.US_state.summary,Province.State %in% top_states_modified ))+
    geom_abline(data=filter(li,Province.State %in% top_states_modified),
                aes(slope = m,intercept = b,col=Province.State),lty=2)+
    geom_point(aes(x=Date,y=Total_confirmed_cases_perstate,col=Province.State))+
    scale_color_brewer(type = "qualitative",palette = "Paired")+
    default_theme+
    theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
    ggtitle("Total confirmed cases per state, top 12 colored")+
    scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Predict the number of total cases over the next week
##-----------------------------------------

predicted_days<-c(0,1,2,3,7)+as.numeric(as.Date("2020-04-20"))

predicted_days_df<-data.frame(matrix(ncol=3))
names(predicted_days_df)<-c("Province.State","days","Total_confirmed_cases_perstate")

# USe model parameters to estiamte case loads
for(state.i in top_states_modified){
  predicted_days_df<-rbind(predicted_days_df,
                           data.frame(Province.State=state.i,
                                      prediction_model(m = li[li$Province.State==state.i,"m"],
                                                       b =li[li$Province.State==state.i,"b"] ,
                                                       days =predicted_days )))
  }

predicted_days_df$Date<-as.Date(predicted_days_df$days,origin="1970-01-01")

kable(predicted_days_df,caption = "Predicted total cases over the next week for selected states")

##------------------------------------------
## Write plots
##-----------------------------------------

write_plot(Corona_Cases.US_state.summary.plot,wd = results_dir)
write_plot(Corona_Cases.US_state.lm.plot,wd = results_dir)

##------------------------------------------
## Write tables
##-----------------------------------------

write.csv(predicted_days_df,file = paste0(results_dir,"predicted_total_cases_days.csv"),quote = F,row.names = F)

Q2: What is the predicted number of cases?

What is the prediction of COVID-19 based on model thus far? Additional questions:

WHy did it take to day 40 to start a log linear trend? How long will it be till x number of cases? When will the plateu happen? Are any effects noticed with social distancing? Delays

##------------------------------------------
## Prediction and Prediction Accuracy
##------------------------------------------


today_num<-max(Corona_Cases.US$Days_since_100)
predicted_days<-today_num+c(1,2,3,7)

#mods = dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)
#today:
Corona_Cases.US[Corona_Cases.US$Days_since_100==(today_num-1),]
Corona_Cases.US[Corona_Cases.US$Days_since_100==today_num,]
Corona_Cases.US$type<-"Historical"


#prediction_values<-prediction_model(m=slope,b=intercept,days = predicted_days)$Total_confirmed_cases

histoical_model<-data.frame(date=today_num,m=slope,b=intercept)
tmp<-data.frame(state=rep(c("A","B"),each=3),x=c(1,2,3,4,5,6))
tmp$y<-c(tmp[1:3,"x"]+5,tmp[4:6,"x"]*5+1)
ddply(tmp,c("state"))
lm(data =tmp,formula = y~x )

train_lm<-function(input_data,subset_coulmn,formula_input){
case_models <- dlply(input_data, subset_coulmn, lm, formula = formula_input)
case_models.parameters <- ldply(case_models, coef)
case_models.parameters<-rename(case_models.parameters,c("b"="(Intercept)","m"=subset_coulmn))
return(case_models.parameters)
}

train_lm(tmp,"state")

 dlply(input_data, subset_coulmn, lm,m=)
 
# model for previous y days
#historical_model_predictions<-data.frame(day_x=NULL,Days_since_100=NULL,Total_confirmed_cases=NULL,Total_confirmed_cases.log=NULL)
# for(i in c(1,2,3,4,5,6,7,8,9,10)){
#   #i<-1
# day_x<-today_num-i # 1, 2, 3, 4
# day_x_nextweek<-day_x+c(1,2,3)
# model_fit_x<-lm(data = filter(Corona_Cases.US.case100,Days_since_100 < day_x),formula = Total_confirmed_cases.log~Days_since_100)
# prediction_day_x_nextweek<-prediction_model(m = model_fit_x$coefficients[2],b = model_fit_x$coefficients[1],days = day_x_nextweek)
# prediction_day_x_nextweek$type<-"Predicted"
# acutal_day_x_nextweek<-filter(Corona_Cases.US,Days_since_100 %in% day_x_nextweek) %>% select(c(Days_since_100,Total_confirmed_cases,Total_confirmed_cases.log))
# acutal_day_x_nextweek$type<-"Historical"
# historical_model_predictions.i<-data.frame(day_x=day_x,rbind(acutal_day_x_nextweek,prediction_day_x_nextweek))
# historical_model_predictions<-rbind(historical_model_predictions.i,historical_model_predictions)
# }

#historical_model_predictions.withHx<-rbind.fill(historical_model_predictions,data.frame(Corona_Cases.US,type="Historical"))
#historical_model_predictions.withHx$Total_confirmed_cases.log2<-log(historical_model_predictions.withHx$Total_confirmed_cases,2)

(historical_model_predictions.plot<-ggplot(historical_model_predictions.withHx,aes(x=Days_since_100,y=Total_confirmed_cases.log,col=type))+
    geom_point(size=3)+
    default_theme+
    theme(legend.position = "bottom")+ 
      #geom_abline(slope = slope,intercept =intercept,lty=2)+
    #facet_wrap(~case_type,ncol=1)+
    scale_color_manual(values = c("Historical"="#377eb8","Predicted"="#e41a1c")))
write_plot(historical_model_predictions.plot,wd=results_dir)

Q3: What is the effect on social distancing, descreased mobility on case load?

Load data from Google which compoutes % change in user mobility relative to baseline for * Recreation
* Workplace
* Residence
* Park
* Grocery

Data from https://www.google.com/covid19/mobility/

# See pre-processing section for script on gathering mobility data

# UNDER DEVELOPMENT

mobility<-read.csv("/Users/stevensmith/Projects/MIT_COVID19/mobility.csv",header = T,stringsAsFactors = F)
#mobility$Retail_Recreation<-as.numeric(sub(mobility$Retail_Recreation,pattern = "%",replacement = ""))
#mobility$Workplace<-as.numeric(sub(mobility$Workplace,pattern = "%",replacement = ""))
#mobility$Residential<-as.numeric(sub(mobility$Residential,pattern = "%",replacement = ""))

##------------------------------------------
## Show relationship between mobility and caseload
##------------------------------------------
mobility$County<-gsub(mobility$County,pattern = " County",replacement = "")
Corona_Cases.US_state.mobility<-merge(Corona_Cases.US_state,plyr::rename(mobility,c("State"="Province.State","County"="City")))

#Corona_Cases.US_state.tmp<-merge(metadata,Corona_Cases.US_state.tmp)
# Needs to happen upsteam, see todos
#Corona_Cases.US_state.tmp$Total_confirmed_cases.perperson<-Corona_Cases.US_state.tmp$Total_confirmed_cases/as.numeric(Corona_Cases.US_state.tmp$Population)
mobility_measures<-c("Retail_Recreation","Grocery_Pharmacy","Parks","Transit","Workplace","Residential")

plot_data<-filter(Corona_Cases.US_state.mobility, Date.numeric==max(Corona_Cases.US_state$Date.numeric) ) %>% melt(measure.vars=mobility_measures) 
plot_data$value<-as.numeric(gsub(plot_data$value,pattern = "%",replacement = ""))
plot_data<-filter(plot_data,!is.na(value))

(mobility.plot<-ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_grid(Province.State~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases per 100 people(Today)"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

(mobility.global.plot<-ggplot(plot_data,aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_wrap(~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases (Today) per 100 people"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

plot_data.permobility_summary<-ddply(plot_data,c("Province.State","variable"),summarise,cor=cor(y =Total_confirmed_cases.per100,x=value),median_change=median(x=value)) %>% arrange(-abs(cor))

kable(plot_data.permobility_summary,caption = "Ranked per-state mobility correlation with total confirmed cases")

Ranked per-state mobility correlation with total confirmed cases

Province.State	variable	cor	median_change
Alaska	Transit	-1.0000000	-63.0
Delaware	Retail_Recreation	1.0000000	-39.5
Delaware	Grocery_Pharmacy	1.0000000	-17.5
Delaware	Parks	-1.0000000	20.5
Delaware	Transit	1.0000000	-37.0
Delaware	Workplace	1.0000000	-37.0
Delaware	Residential	-1.0000000	14.0
Hawaii	Grocery_Pharmacy	0.9996038	-34.0
Hawaii	Retail_Recreation	0.9947437	-56.0
New Hampshire	Parks	0.9533349	-20.0
Connecticut	Grocery_Pharmacy	-0.9046327	-6.0
Alaska	Residential	0.8967338	13.0
South Dakota	Parks	0.8909933	-26.0
Maine	Transit	-0.8885916	-50.0
Vermont	Parks	0.8581722	-35.5
Wyoming	Transit	-0.8418991	-16.0
Hawaii	Parks	0.8242552	-72.0
Utah	Residential	-0.8170563	12.0
Alaska	Grocery_Pharmacy	-0.8134720	-7.0
Utah	Workplace	-0.8129936	-37.0
Massachusetts	Workplace	-0.7851536	-39.0
Utah	Parks	-0.7794775	17.0
Connecticut	Transit	-0.7781216	-50.0
Hawaii	Transit	0.7749229	-89.0
Utah	Transit	-0.7595546	-18.0
Rhode Island	Workplace	-0.7474680	-39.5
Alaska	Workplace	-0.7214938	-34.0
Maine	Workplace	-0.6723722	-30.0
Vermont	Grocery_Pharmacy	-0.6520073	-25.0
New York	Workplace	-0.6417529	-34.5
New Jersey	Workplace	-0.6391067	-44.0
Montana	Workplace	-0.6351999	-40.5
Hawaii	Residential	-0.6320557	19.0
Rhode Island	Residential	-0.6263357	18.5
Rhode Island	Retail_Recreation	-0.6088746	-45.0
Nebraska	Workplace	0.6081270	-32.5
Arizona	Grocery_Pharmacy	-0.5908662	-15.0
New York	Retail_Recreation	-0.5836636	-46.0
North Dakota	Retail_Recreation	-0.5795886	-43.5
Wyoming	Parks	-0.5791942	-4.0
New Jersey	Parks	-0.5775734	-6.0
Connecticut	Residential	0.5599949	14.0
Massachusetts	Retail_Recreation	-0.5505991	-44.0
West Virginia	Parks	0.5413903	-33.0
New Jersey	Retail_Recreation	-0.5407021	-62.5
Maine	Parks	0.5284362	-31.0
New York	Parks	0.5208609	20.0
Arkansas	Parks	-0.5085784	-12.0
Connecticut	Workplace	-0.5013570	-39.0
Nebraska	Residential	-0.4956025	14.0
Montana	Parks	-0.4913929	-58.0
Wyoming	Workplace	-0.4894702	-31.0
Iowa	Parks	-0.4889626	28.5
New Jersey	Grocery_Pharmacy	-0.4877835	2.5
Connecticut	Retail_Recreation	-0.4833003	-45.0
Utah	Retail_Recreation	-0.4825461	-40.0
New Mexico	Grocery_Pharmacy	-0.4803291	-11.0
Arizona	Retail_Recreation	-0.4788668	-42.5
Rhode Island	Parks	0.4782766	52.0
Idaho	Workplace	-0.4740073	-29.0
Wisconsin	Transit	-0.4729892	-23.5
Montana	Residential	0.4679731	14.0
New Mexico	Parks	0.4635243	-31.5
Illinois	Transit	-0.4608843	-31.0
Massachusetts	Grocery_Pharmacy	-0.4521142	-7.0
Idaho	Transit	-0.4480548	-30.0
New Mexico	Residential	0.4464201	13.5
California	Transit	-0.4457163	-42.0
California	Residential	0.4425696	14.0
Montana	Retail_Recreation	-0.4329731	-51.0
New Jersey	Transit	-0.4300372	-50.5
Pennsylvania	Workplace	-0.4298672	-36.0
Vermont	Residential	0.4281670	11.5
Kentucky	Parks	-0.4252020	28.5
New Hampshire	Residential	-0.4177960	14.0
South Carolina	Workplace	0.4171950	-30.0
Idaho	Grocery_Pharmacy	-0.3971667	-4.5
Montana	Grocery_Pharmacy	-0.3906632	-16.0
Montana	Transit	-0.3892824	-41.0
Alabama	Workplace	-0.3855478	-29.0
Alabama	Transit	-0.3817654	-36.5
Maryland	Grocery_Pharmacy	-0.3759528	-10.0
New Mexico	Retail_Recreation	-0.3747372	-42.0
New York	Transit	-0.3741999	-48.0
Florida	Residential	0.3698720	14.0
North Dakota	Workplace	0.3655293	-38.0
California	Parks	-0.3604178	-38.5
Maryland	Workplace	-0.3590174	-35.0
Pennsylvania	Retail_Recreation	-0.3540238	-45.0
Wyoming	Grocery_Pharmacy	-0.3437261	-9.0
Hawaii	Workplace	0.3435834	-46.0
Kansas	Parks	0.3420112	72.0
Arizona	Residential	0.3418058	13.0
Idaho	Retail_Recreation	-0.3387792	-40.5
Nebraska	Grocery_Pharmacy	0.3318352	-0.5
Arizona	Transit	0.3314223	-38.0
Alabama	Grocery_Pharmacy	-0.3288560	-2.0
California	Retail_Recreation	-0.3287409	-44.0
Nevada	Transit	-0.3286630	-20.0
Pennsylvania	Parks	0.3202615	13.0
Vermont	Retail_Recreation	0.3200580	-57.0
North Dakota	Parks	0.3193606	-34.0
Minnesota	Transit	-0.3184238	-28.5
Alaska	Retail_Recreation	0.3155306	-39.0
California	Workplace	-0.3130160	-36.0
Florida	Parks	-0.3073586	-43.0
California	Grocery_Pharmacy	-0.3042905	-12.0
Michigan	Parks	0.3024387	30.0
Colorado	Residential	0.2982059	14.0
Maine	Retail_Recreation	-0.2934388	-42.0
Arkansas	Retail_Recreation	-0.2903966	-30.0
North Carolina	Grocery_Pharmacy	0.2879048	0.0
Kansas	Workplace	0.2851573	-33.0
Texas	Workplace	0.2831242	-32.0
Rhode Island	Grocery_Pharmacy	0.2822096	-7.5
Rhode Island	Transit	-0.2807816	-56.0
Nevada	Retail_Recreation	-0.2803360	-43.0
Texas	Residential	-0.2779958	15.0
Virginia	Transit	-0.2747326	-33.0
Oregon	Grocery_Pharmacy	0.2745336	-7.0
Nevada	Residential	0.2696828	17.0
Tennessee	Workplace	-0.2684606	-31.0
Tennessee	Residential	0.2667282	11.5
Maryland	Retail_Recreation	-0.2659919	-39.0
Texas	Parks	0.2639738	-42.0
West Virginia	Grocery_Pharmacy	-0.2566666	-6.0
Wisconsin	Parks	0.2530271	51.5
Georgia	Grocery_Pharmacy	-0.2528004	-10.0
Illinois	Workplace	-0.2525636	-31.0
North Dakota	Grocery_Pharmacy	-0.2524942	-8.0
Maryland	Residential	0.2442600	15.0
North Carolina	Workplace	0.2429153	-31.5
Michigan	Workplace	-0.2421589	-40.0
Mississippi	Residential	0.2412199	13.0
North Carolina	Transit	0.2395306	-32.0
Alabama	Parks	0.2390554	-1.0
New York	Grocery_Pharmacy	-0.2336830	8.0
Pennsylvania	Grocery_Pharmacy	-0.2307583	-6.0
North Carolina	Residential	0.2260297	13.0
Arkansas	Residential	0.2199946	12.0
Iowa	Transit	0.2192411	-24.0
South Dakota	Transit	-0.2161209	-40.0
New Jersey	Residential	0.2140058	18.0
Maine	Grocery_Pharmacy	-0.2111839	-13.0
West Virginia	Workplace	0.2109873	-33.0
South Carolina	Parks	-0.2091405	-23.0
Oregon	Residential	0.2084809	10.5
Washington	Workplace	-0.2082184	-38.0
Utah	Grocery_Pharmacy	-0.2075200	-4.0
Kansas	Grocery_Pharmacy	-0.2040020	-14.0
Illinois	Parks	0.2018857	26.5
Mississippi	Grocery_Pharmacy	-0.2005005	-8.0
Georgia	Workplace	-0.1996536	-33.5
Wyoming	Retail_Recreation	-0.1932005	-40.0
Texas	Transit	0.1905799	-42.0
Colorado	Parks	-0.1903552	2.0
Missouri	Workplace	0.1856128	-28.5
Virginia	Grocery_Pharmacy	-0.1845875	-8.0
Oklahoma	Parks	-0.1806422	-18.5
Vermont	Workplace	-0.1794875	-43.0
Georgia	Retail_Recreation	-0.1788391	-41.0
Washington	Transit	-0.1762795	-33.5
Massachusetts	Transit	-0.1761902	-45.0
Illinois	Residential	0.1720879	14.0
Nebraska	Parks	0.1698799	55.5
Maine	Residential	-0.1694094	11.0
Michigan	Retail_Recreation	-0.1684882	-53.0
Florida	Workplace	-0.1671606	-33.0
Ohio	Transit	0.1670829	-28.0
Virginia	Parks	0.1669170	6.0
Oklahoma	Residential	0.1667877	15.0
New Mexico	Transit	0.1661050	-38.0
Missouri	Residential	-0.1643004	13.0
Virginia	Residential	0.1608895	14.0
Georgia	Residential	-0.1590277	13.0
Wisconsin	Residential	-0.1578373	14.0
South Dakota	Retail_Recreation	-0.1528491	-38.5
New Hampshire	Transit	-0.1511516	-57.0
Massachusetts	Residential	0.1475725	15.0
West Virginia	Residential	-0.1470979	11.0
South Carolina	Residential	-0.1454668	12.0
Connecticut	Parks	0.1420883	43.0
Oregon	Parks	0.1409019	16.5
Ohio	Parks	-0.1392998	67.5
Michigan	Grocery_Pharmacy	-0.1350377	-11.0
Arkansas	Workplace	-0.1344201	-26.0
Indiana	Residential	0.1334621	12.0
Indiana	Retail_Recreation	0.1318773	-38.0
Idaho	Parks	0.1317838	-22.0
Washington	Residential	0.1311443	13.0
West Virginia	Retail_Recreation	0.1306241	-38.5
North Carolina	Parks	-0.1288207	7.0
Mississippi	Workplace	-0.1268884	-33.0
Mississippi	Retail_Recreation	-0.1263908	-40.0
Pennsylvania	Transit	-0.1261784	-41.5
Florida	Retail_Recreation	0.1256920	-43.0
Minnesota	Workplace	-0.1247611	-33.0
Oregon	Retail_Recreation	0.1225698	-41.0
Massachusetts	Parks	0.1213375	39.0
New Hampshire	Retail_Recreation	-0.1213289	-41.0
Kansas	Transit	-0.1212614	-26.5
Arkansas	Transit	0.1200796	-27.0
Maryland	Transit	-0.1151047	-39.0
Nebraska	Retail_Recreation	0.1149303	-36.0
Kentucky	Transit	0.1137272	-31.0
Mississippi	Transit	-0.1132114	-38.5
Indiana	Parks	-0.1129299	29.0
Washington	Parks	0.1125132	-3.5
Alabama	Retail_Recreation	0.1110507	-39.0
New Hampshire	Grocery_Pharmacy	-0.1104435	-6.0
Washington	Grocery_Pharmacy	0.1079446	-7.0
Ohio	Residential	0.1067426	14.0
Pennsylvania	Residential	0.1028482	15.0
Wisconsin	Grocery_Pharmacy	0.1020410	-1.0
Kentucky	Grocery_Pharmacy	0.1015358	4.0
South Dakota	Residential	0.1008153	15.0
Missouri	Transit	-0.0995212	-24.5
Arizona	Workplace	-0.0978333	-35.0
South Dakota	Grocery_Pharmacy	0.0974784	-9.0
Idaho	Residential	-0.0967765	11.0
Oregon	Workplace	-0.0945299	-31.0
Indiana	Workplace	0.0937960	-34.0
Wyoming	Residential	0.0910868	12.5
New York	Residential	0.0906452	17.5
Wisconsin	Workplace	-0.0887610	-31.0
Iowa	Retail_Recreation	-0.0876596	-38.0
Virginia	Workplace	-0.0842931	-31.5
South Carolina	Transit	0.0826127	-45.0
Georgia	Parks	0.0811019	-6.0
Minnesota	Retail_Recreation	0.0797620	-40.0
Virginia	Retail_Recreation	-0.0795494	-35.0
Tennessee	Transit	-0.0791873	-32.0
Minnesota	Parks	0.0775662	-9.0
South Carolina	Retail_Recreation	-0.0761711	-35.0
Texas	Grocery_Pharmacy	0.0754527	-14.0
Indiana	Grocery_Pharmacy	-0.0750331	-5.5
Oklahoma	Grocery_Pharmacy	-0.0744429	-1.0
Ohio	Retail_Recreation	0.0737054	-36.0
Michigan	Residential	0.0733463	15.0
Kentucky	Retail_Recreation	0.0701951	-29.0
Iowa	Workplace	-0.0689782	-30.5
Ohio	Grocery_Pharmacy	0.0686977	0.0
Florida	Transit	-0.0664891	-49.0
Kentucky	Residential	0.0636959	12.0
Colorado	Transit	0.0629383	-36.0
Colorado	Grocery_Pharmacy	-0.0604831	-17.0
Colorado	Retail_Recreation	-0.0596211	-44.0
New Hampshire	Workplace	0.0586084	-37.0
Oregon	Transit	0.0543691	-27.5
Illinois	Retail_Recreation	0.0525875	-40.0
Minnesota	Grocery_Pharmacy	0.0509333	-6.5
Washington	Retail_Recreation	-0.0483603	-42.0
Kentucky	Workplace	-0.0477633	-36.0
Michigan	Transit	0.0473397	-46.0
Tennessee	Parks	-0.0472954	10.5
Nevada	Workplace	0.0472344	-40.0
Texas	Retail_Recreation	0.0459951	-40.0
South Carolina	Grocery_Pharmacy	0.0432394	1.0
North Dakota	Transit	0.0428862	-48.0
Missouri	Retail_Recreation	-0.0402379	-36.0
Missouri	Parks	0.0372349	0.0
North Dakota	Residential	0.0341653	17.0
Missouri	Grocery_Pharmacy	0.0338297	2.0
Oklahoma	Retail_Recreation	0.0332175	-31.0
West Virginia	Transit	0.0314701	-45.0
Minnesota	Residential	0.0296641	17.0
Illinois	Grocery_Pharmacy	-0.0291781	2.0
Tennessee	Grocery_Pharmacy	0.0271857	6.0
Colorado	Workplace	-0.0259667	-39.0
Kansas	Residential	-0.0257811	13.0
Oklahoma	Workplace	0.0256221	-31.0
Vermont	Transit	0.0251328	-63.0
Alabama	Residential	0.0239661	11.0
Arkansas	Grocery_Pharmacy	-0.0214689	3.0
Nebraska	Transit	-0.0201795	-9.0
Iowa	Grocery_Pharmacy	0.0199525	4.0
Florida	Grocery_Pharmacy	0.0191262	-14.0
North Carolina	Retail_Recreation	0.0186505	-34.0
Wisconsin	Retail_Recreation	0.0174096	-44.5
Kansas	Retail_Recreation	-0.0158722	-38.0
Ohio	Workplace	-0.0136880	-35.0
Arizona	Parks	-0.0125018	-44.5
Nevada	Parks	-0.0112991	-12.5
Maryland	Parks	-0.0088106	27.0
Tennessee	Retail_Recreation	-0.0059371	-30.0
Iowa	Residential	-0.0044421	13.0
New Mexico	Workplace	0.0038019	-34.0
Oklahoma	Transit	0.0030317	-26.0
Nevada	Grocery_Pharmacy	0.0029789	-12.5
Georgia	Transit	-0.0023010	-35.0
Mississippi	Parks	0.0021238	-25.0
South Dakota	Workplace	0.0015186	-35.0
Indiana	Transit	-0.0003852	-29.0
Alaska	Parks	NA	29.0
District of Columbia	Retail_Recreation	NA	-69.0
District of Columbia	Grocery_Pharmacy	NA	-28.0
District of Columbia	Parks	NA	-65.0
District of Columbia	Transit	NA	-69.0
District of Columbia	Workplace	NA	-48.0
District of Columbia	Residential	NA	17.0

# sanity check
ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(x=Total_confirmed_cases.per100,fill=variable))+geom_histogram()+
  facet_grid(~Province.State)+
    default_theme+
  theme(legend.position = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

write_plot(mobility.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.plot.png"

write_plot(mobility.global.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.global.plot.png"

(plot_data.permobility_summary.plot<-ggplot(plot_data.permobility_summary,aes(x=variable,y=median_change))+
  geom_jitter(size=2,width=.2)+
  #geom_jitter(data=plot_data.permobility_summary %>% arrange(-abs(median_change)) %>% head(n=15),aes(col=Province.State),size=2,width=.2)+
  default_theme+
  ggtitle("Per-Sate Median Change in Mobility")+
  xlab("Mobility Meaure")+
  ylab("Median Change from Baseline"))

write_plot(plot_data.permobility_summary.plot,wd = results_dir)

## [1] "/Users/stevensmith/Projects/coronavirus/results/plot_data.permobility_summary.plot.png"

DELIVERABLE MANIFEST

The following link to commited documents pushed to github. These are provided as a convienence, but note this is a manual process. The generation of reports, plots and tables is not coupled to the execution of this markdown. ## Report This report, html & pdf

Plots

github_root<-"https://github.com/sbs87/coronavirus/blob/master/"

plot_handle<-c("Corona_Cases.world.long.plot",
               "Corona_Cases.world.loglong.plot",
               "Corona_Cases.world.mortality.plot",
               "Corona_Cases.world.casecor.plot",
               "Corona_Cases.city.long.plot",
               "Corona_Cases.city.loglong.plot",
               "Corona_Cases.city.mortality.plot",
               "Corona_Cases.city.casecor.plot",
               "Corona_Cases.city.long.normalized.plot",
               "Corona_Cases.US_state.lm.plot",
               "Corona_Cases.US_state.summary.plot")


deliverable_manifest<-data.frame(
  name=c("World total & death cases, longitudinal",
         "World log total & death cases, longitudinal",
         "World mortality",
         "World total & death cases, correlation",
         "City total & death cases, longitudinal",
         "City log total & death cases, longitudinal",
         "City mortality",
         "City total & death cases, correlation",
         "City population normalized total & death cases, longitudinal",
         "State total cases (select) with linear model, longitudinal",
         "State total cases, longitudinal"),
  plot_handle=plot_handle,
  link=paste0(github_root,"results/",plot_handle,".png")
)


(tmp<-data.frame(row_out=apply(deliverable_manifest,MARGIN = 1,FUN = function(x) paste(x[1],x[2],x[3],sep=" | "))))

##                                                                                                                                                                                                        row_out
## 1                                           World total & death cases, longitudinal | Corona_Cases.world.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
## 2                                 World log total & death cases, longitudinal | Corona_Cases.world.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
## 3                                                         World mortality | Corona_Cases.world.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
## 4                                      World total & death cases, correlation | Corona_Cases.world.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
## 5                                              City total & death cases, longitudinal | Corona_Cases.city.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
## 6                                    City log total & death cases, longitudinal | Corona_Cases.city.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
## 7                                                            City mortality | Corona_Cases.city.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
## 8                                         City total & death cases, correlation | Corona_Cases.city.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
## 9  City population normalized total & death cases, longitudinal | Corona_Cases.city.long.normalized.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
## 10                     State total cases (select) with linear model, longitudinal | Corona_Cases.US_state.lm.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
## 11                                      State total cases, longitudinal | Corona_Cases.US_state.summary.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

row_out<-apply(tmp, 2, paste, collapse="\t\n")

name	handle	link
World total & death cases, longitudinal	Corona_Cases.world.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
World log total & death cases, longitudinal	Corona_Cases.world.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
World mortality	Corona_Cases.world.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
World total & death cases, correlation	Corona_Cases.world.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
City total & death cases, longitudinal	Corona_Cases.city.long.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
City log total & death cases, longitudinal	Corona_Cases.city.loglong.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
City mortality	Corona_Cases.city.mortality.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
City total & death cases, correlation	Corona_Cases.city.casecor.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
City population normalized total & death cases, longitudinal	Corona_Cases.city.long.normalized.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
State total cases (select) with linear model, longitudinal	Corona_Cases.US_state.lm.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
State total cases, longitudinal	Corona_Cases.US_state.summary.plot	https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

Tables

CONCLUSION

Overall, the trends of COVID-19 cases is no longer in log-linear phase for world or U.S. (but some regions like MD are still in the log-linear phase). Mortality rate (deaths/confirmed RNA-based cases) is >1%, with a range depending on region. Mobility is not a strong indicator of caseload (U.S. data).

See table below for detailed breakdown.

Question	Answer
What is the effect on social distancing, descreased mobility on case load?	There is not a strong apparent effect on decreased mobility (work, grocery, retail) or increased mobility (at residence, parks) on number of confirmed cases, either as a country (U.S.) or state level. California appears to have one of the best correlations, but this is a mixed bag
What is the trend in cases, mortality across geopgraphical regions?	The confirmed total casees and mortality is overall log-linear for most countries, with a trailing off beginning for most (inlcuding U.S.). On the state level, NY, NJ, PA starting to trail off; MD is still in log-linear phase. Mortality and case load are highly correlated for NY, NJ, PA, MD. The mortality rate flucutates for a given region, but is about 3% overall.

END

End: ##—— Sun May 17 19:32:01 2020 ——##

Cheatsheet: http://rmarkdown.rstudio.com>

Sandbox

# Geographical heatmap!
install.packages("maps")
library(maps)
library
mi_counties <- map_data("county", "pennsylvania") %>% 
  select(lon = long, lat, group, id = subregion)
head(mi_counties)

ggplot(mi_counties, aes(lon, lat)) + 
  geom_point(size = .25, show.legend = FALSE) +
  coord_quickmap()
mi_counties$cases<-1:2226
name_overlaps(metadata,Corona_Cases.US_state)

tmp<-merge(Corona_Cases.US_state,metadata)
ggplot(filter(tmp,Province.State=="Pennsylvania"), aes(Long, Lat, group = as.factor(City))) +
  geom_polygon(aes(fill = Total_confirmed_cases), colour = "grey50") + 
  coord_quickmap()


ggplot(Corona_Cases.US_state, aes(Long, Lat))+
  geom_polygon(aes(fill = Total_confirmed_cases ), color = "white")+
  scale_fill_viridis_c(option = "C")
dev.off()


require(maps)
require(viridis)

world_map <- map_data("world")
ggplot(world_map, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill="lightgray", colour = "white")

head(world_map)
head(Corona_Cases.US_state)
unique(select(world_map,c("region","group"))) %>% filter()

some.eu.countries <- c(
  "US"
)
# Retrievethe map data
some.eu.maps <- map_data("world", region = some.eu.countries)

# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region.lab.data <- some.eu.maps %>%
  group_by(region) %>%
  summarise(long = mean(long), lat = mean(lat))

unique(filter(some.eu.maps,subregion %in% Corona_Cases.US_state$Province.State) %>% select(subregion))
unique(Corona_Cases.US_state$Total_confirmed_cases.log)
ggplot(filter(Corona_Cases.US_state,Date=="2020-04-17") aes(x = Long, y = Lat)) +
  geom_polygon(aes( fill = Total_confirmed_cases.log))+
  #geom_text(aes(label = region), data = region.lab.data,  size = 3, hjust = 0.5)+
  #scale_fill_viridis_d()+
  #theme_void()+
  theme(legend.position = "none")
library("sf")
library("rnaturalearth")
library("rnaturalearthdata")

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
ggplot(data = world) +
    geom_sf()

counties <- st_as_sf(map("county", plot = FALSE, fill = TRUE))
counties <- subset(counties, grepl("florida", counties$ID))
counties$area <- as.numeric(st_area(counties))
#install.packages("lwgeom")
class(counties)
head(counties)
ggplot(data = world) +
    geom_sf(data=Corona_Cases.US_state) +
    #geom_sf(data = counties, aes(fill = area)) +
  geom_sf(data = counties, aes(fill = area)) +
   # scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
    coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)


head(counties)
tmp<-unique(select(filter(Corona_Cases.US_state,Date=="2020-04-17"),c(Lat,Long,Total_confirmed_cases.per100)))
st_as_sf(map("county", plot = FALSE, fill = TRUE))

join::inner_join.sf(Corona_Cases.US_state, counties)

library(sf)
library(sp)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
class(nc)


spdf <- SpatialPointsDataFrame(coords = select(Corona_Cases.US_state,c("Lat","Long")), data = Corona_Cases.US_state,
                               proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

head(spdf)
class(spdf)
st_cast(spdf)

filter(Corona_Cases.US_state.summary,Date=="2020-04-20" & Province.State %in% top_states_modified)

https://stevenbsmith.net